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I.   INTRODUCTION 

A  most  critical  portion  of  flight,  even  in  an  ideal 
operating  environment,  is  the  landing  of  an  airplane  at  a 
desired  landing  site.   Inclement  weather  at  the  landing 
site,  with  the  corresponding  reduction  in  visibility,  fur- 
ther complicates  this  portion  of  the  flight  operation  to 
the  point  that  landing  may  become  an  impossibility  due  to 
visibility  restrictions.   The  increased  demand,  both  com- 
mercially and  militarily,  for  flight  operations  under  all- 
weather  conditions  has  prompted  a  great  deal  of  research 
to  develop  an  "Automatic  Ail-Weather  Landing  System"  capa- 
ble of  automatically  controlling  an  airplane  during  the 
landing  portion  of  a  flight.   The  problems  encountered  in 
such  investigations,  although  not  considered  insurmount- 
able, are  vast  and  complicated  mainly  because  the  flight 
dynamics  of  the  airplane  and  the  surrounding  environment 
are  themselves  complicated.   This  investigation  makes  some 
simplifying  assumptions  to  reduce  this  complexity  and  de- 
rives the  necessary  controls  for  such  an  all-weather  land- 
ing system,  using  optimal  control  techniques. 

The  problem  and  the  simplifying  assumptions  are  de- 
scribed in  Chapter  II.   In  Chapter  III  a  general  model  for 
an  airplane  in  the  landing  mode  is  developed,  and  the  spe- 
cific model  used  in  this  investigation  is  presented. 
Chapter  IV  presents  the  optimal  control  theory  employed 
and  describes  the  method  used  to  modify  the  original  state 


variables  to  make  the  theory  applicable.   The  problem 
specifications  are  introduced  in  Chapter  V,  the  desired 
trajectory  is  formulated,  and  the  significance  of  the  per- 
formance measure  in  obtaining  realistic  results  is  dis- 
cussed.  In  Chapter  VI  the  procedures  used  in  the  investi- 
gation are  explained,  and  the  cases  investigated  are  pre- 
sented.  The  results  are  discussed  in  Chapter  VII,  and  the 
conclusions  and  recommendations  are  presented  in  Chapter 
VIII. 
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II.   PROBLEM  FORMULATION 

The  basic  descriptions  and  definitions  of  the  landing 
portion  of  flight  are  presented  in  this  chapter,  along  with 
the  assumptions  made  to  reduce  the  complexity  of  the  pro- 
blem. 

A.    STANDARD  WEATHER  CRITERIA 

The  International  Civil  Aviation  Organization  (ICAO) 
has  adopted  the  following  set  of  minimum  weather  conditions 
for  the  automatic  landing  of  an  airplane,  (Ref .  1) .   All 
distance  specifications  are  given  in  meters,  followed  in 
parenthesis  by  approximate  values  in  feet. 

Category  I.   Operation  down  to  minima  of  60  meters 
(200  ft.)  decision  height  (altitude)  and  Runway  Visibility 
Reading  (RVR)  of  800  meters  (2600  ft.). 

Category  II.   Operation  down  to  minima  of  30  meters 
(100  ft.)  decision  height  and  a  RVR  of  400  meters  (1200  ft.). 

Category  III-A.   Operation  to  and  along  the  surface  of 
the  runway,  with  external  visual  reference  during  the  final 
phase  of  the  landing  to  a  RVR  minimum  of  200  meters  (700 
ft.)  . 

Category  III-B.   Operation  to  and  along  the  surface  of 
the  runway  and  taxiways,  with  visibility  sufficient  only 
for  visual  taxiing  comparable  to  an  RVR  of  about  50  meters 
(150  ft.)  . 
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Category  III-C.   Operation  to  and  along  the  surface  of 
the  runway  and  taxiway  without  external  visual  reference. 

The  Federal  Aviation  Agency  (FAA)  has  established  the 
criteria  for  certification  of  any  automatic  all-weather 
landing  system  based  on  compliance  with  these  categories  of 
automatic  landing  weather  minima.   Presently  no  automatic 
systems  that  comply  with  any  Category  III  minima  have  been 
certified.   However,  experimental  tests  have  been  success- 
fully conducted  (Category  III-B)  with  a  system  in  a  C-141 
airplane  (Ref.  2).   Several  automatic  control  systems  have 
been  certified  for  Category  II  operation,  provided  that  the 
landing  site  has  the  prerequisite  ground  equipment.   The 
vast  majority  of  automatic  landings  are  restricted  to 
Category  I  conditions. 

B.    GROUND  EQUIPMENT  AT  AN  ALL-WEATHER  LANDING  SITE 

At  present,  only  one  special  piece  of  ground  equipment, 
an  Instrument  Landing  Systems  (ILS) ,  is  required  for  land- 
ing at  an  all-weather  landing  facility.   The  ILS  provides  a 
radio  beam  to  establish  a  nominal  glide  path  to  the  land- 
ing site  and  several  locator  beacons,  placed  along  the  land- 
ing track  to  provide  distance  checks.   One  locator  beacon 
is  positioned  to  provide  an  indication  of  the  nominal  deci- 
sion altitude  point  for  the  facility.   Both  azimuth  and 
elevation  information  are  provided  by  the  ILS.   Category  I 
and  Category  II  automatic  landing  operations  are  presently 
using  this  ground  equipment  with  the  latter,  of  necessity, 
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requiring  more  sophisticated  ILS  equipment  to  provide  the 
increased  accuracy  needed  in  the  radio  beam.   Figure  1 
shows  present  requirements  for  the  beam  accuracy  of  ILS 
Category-II  ground  equipment.   Several  major  U.S.  air- 
fields have  this  equipment  installed  and  operational. 
Because  of  beam  distortions  due  to  the  electromagnetic  dis- 
turbances, antenna  installation  environments,  etc.,  the  ILS 
equipment  is  not  accurate  below  an  altitude  of  about  60 
feet.   Therefore,  any  automatic  control  of  the  airplane 
beyond  Category  II  minima  will  require  either  additional 
ground  equipment  or  an  automatic  system  in  the  airplane  it- 
self to  effect  the  final  phases  of  the  landing. 


1100  f 


ii__i_ 


Figure  1 
Category  II  ILS  Beam  Accuracy  Requirements 
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C.    PHASES  OF  ALL-WEATHER  LANDING 

The  automatic  landing  of  an  airplane  under  all-weather 
conditions  can  be  divided  into  five  phases: 

1.  The  Approach-Intercept  Phase 

The  Approach-Intercept  Phase  is  defined  as  that 
segment  of  the  landing  during  which  the  configuration  of 
the  airplane  is  physically  changed  in  anticipation  of  the 
landing  and  the  airplane  is  directed  toward  the  ILS  beam 
which  serves  the  landing  site. 

2.  The  Capture-Track  Phase 

The  Capture-Track  Phase  is  defined  as  that  segment  of 
the  landing  during  which  the  airplane  is  directed  toward 
the  landing  point  by  an  automatic  flight  control  system 
which  uses  the  ILS  beam  as  its  reference. 

3 .  The  Flare  Phase 

The  Flare  Phase  is  defined  as  that  segment  of  the 
landing  during  which  the  airplane  attitude  is  changed  in 
preparation  for  the  touchdown. 

4.  The  Land  Phase 

The  Land  Phase  is  defined  as  the  actual  touchdown 
of  the  airplane  on  the  landing  surface. 

5.  The  Roll-Out  Phase 

The  Roll-Out  Phase  is  defined  as  the  climax  of 
the  landing  during  which  the  airplane  is  slowed  to  a  stop 
or,  equivalently,  to  a  safe  taxi  speed. 
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Figure  2  depicts  the  phases  of  landing  in  both  azimuth 
and  pitch.   Included,  for  clarity,  is  the  decision  altitude 
(height) ,  which  is  defined  as  the  altitude  at  which  a  deci- 
sion must  be  made  as  to  whether  it  is  feasible  to  continue 
the  landing.   As  an  example,  the  decision  altitude  for 
Category  II  landings  is  100  feet  above  ground  level. 

Automatic  control  of  the  airplane  during  the  first 
three  phases  of  the  landing  is  presently  standard  operating 
procedure  for  both  commercial  and  military  flight  opera- 
tions.  However,  lack  of  certified  automatic  flight  control 
systems  to  accomplish  the  remaining  phases  of  the  landing 
dictates  that  the  decision  to  continue  the  landing  beyond 
the  decision  altitude  be  based  on  the  premise  that  the 
pilot  has  visual  contact  with  the  landing  site  at  the  de- 
cision altitude,  so  that  the  airplane  can  be  controlled 
manually,  if  necessary. 

D.    ASSUMPTIONS  MADE  FOR  THE  LANDING  PROBLEM  IN  THIS  STUDY 
As  stated  in  the  opening  remarks  of  this  section,  some 
simplifying  assumptions  are  made  in  this  study  to  reduce  the 
complexity  of  the  automatic  landing  problem.   These  assump- 
tions are  listed  below. 

Assumption  1.   Only  that  portion  of  the  landing  from 
the  decision  altitude  to  the  land  phase  will  be  considered. 
It  is  assumed  that  the  airplane  is  automatically  or  manually 
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controlled,  using  the  ILS  as  a  reference,  to  the  locator 
beacon  which  defines  the  nominal  decision  altitude  point  -- 
taken  to  be  100  feet  above  ground  level.   The  roll-out 
phase  is  also  assumed  to  be  automatically  or  manually  con- 
trolled. 

Assumption  2.   The  airplane  is  physically  located  in 
space  within  the  prescribed  Category-II  ILS  "window"  at  the 
nominal  decision  altitude  point  and  is  in  Equilibrium 
Flight  (see  definition  in  Chapter  III)  at  that  time.   If 
these  conditions  are  not  met,  it  is  assumed  that  the  land- 
ing will  be  discontinued  (airplane  will  be  waved  off) . 

Assumption  3.   Wind  effects  will  not  be  considered. 
During  landing,  the  airplane  is  subjected  to  both  steady- 
state  and  gusty  winds,  the  latter  being  of  primary  impor- 
tance, since  the  gusts  are  random  in  nature   Steady-state 
winds  could  be  considered  since  their  effect  can  be  elimi- 
nated by  a  steady-state  change  in  the  airplane  heading, 
whereas  the  random  wind  gusts  could  require  a  statistical 
handling  of  the  problem.   By  neglecting  both  effects,  the 
airplane  velocity  becomes  equal  to  the  ground  velocity  and 
the  problem  can  be  formulated  in  terms  of  time-to-go-to 
landing. 

Assumption  4.   Only  airplane  motion  in  the  vertical 
plane  will  be  considered.   Lateral  motion  of  the  airplane 
in  the  final  phases  of  landing  is  primarily  necessary  to 
point  the  airplane  in  the  direction  of  the  runway  just  prior 
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to  the  actual  touchdown.   The  airplane  is  normally  "off" 
heading  during  the  landing  to  counter  steady-state  cross 
winds.   Since  the  wind  effects  have  been  neglected,  the 
lateral  motion  can  be  neglected. 
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III.   MODELING  THE  PLANT 

The  first  step  in  any  control  problem  is  the  formula- 
tion of  a  realistic  mathematical  model  to  represent  the 
dynamics  of  the  plant  to  be  controlled  --  in  this  case,  the 
airplane.   In  this  chapter  the  formualtion  of  a  general 
model  applicable  to  any  airplane  is  discussed,  and  the  spe- 
cific model  used  in  this  investigation  is  presented. 

A.    THE  GENERAL  MATHEMATICAL  MODEL 

Aerodynamists  have  developed  two  sets  of  linear  equa- 
tions which  describe  the  dynamics  of  an  airplane.   These 
are  referred  to  as  the  longitudinal  or  symmetric  and  the 
lateral  or  asymmetric  equations  of  motion.   Since,  as  pre- 
viously stated,  this  problem  considers  only  the  longitudi- 
nal motions  of  the  airplane,  further  reference  is  limited 
to  the  longitudinal  equations  of  motion  only.   Derivation 
of  these  equations  is  beyond  the  scope  of  this  study,  but 
complete  and  detailed  derivations  can  be  found  in  the  li- 
terature.  (Ref.  3,  4,  and  5) 

Before  presenting  the  equations  of  motion,  it  is  sig- 
nificant, for  clarification  and  reference  purposes,  to  pre- 
sent several  aerodynamic  definitions  and  assumptions  used 
in  obtaining  these  equations. 

1.   Aerodynamic  Definitions 

a.   Equilibrium  flight  is  defined  as  unaccelerated 
flight.   Flight  is  along  a  straight  path  during  which  the 
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linear  velocity  vector  measured  relative  to  fixed  space 
is  invariant  and  angular  velocity  is  zero. 

b.  Steady  flight  is  defined  as  flight  during 
which  the  linear  velocity  vector  is  invariant  and  angular 
velocity  is  constant.   In  this  context,  equilibrium  flight 
is  also  steady  flight. 

c.  Airplane  coordinates,  axes,  and  angles  are 
defined  in  Figure  3,  which  shows  the  airplane  in  equili- 
brium flight. 

v   =  Equilibrium  flight  linear  velocity  along  x  axis 

Oi     =  Equilibrium  flight  angle  of  attack 
0   =  Equilibrium  flight  pitch  angle 
y      ~   Equilibrium  flight  glide  angle 

d.  Instantaneous  components  are  defined  as  the 
summation  of  equilibrium  flight  components  and  correspond- 
ing perturbations  caused  by  a  disturbed  flight  condition. 
Figure  4  is  a  representation  of  the  airplane  in  a  disturbed 
flight  condition  and  includes  the  perturbations  of  the  two 
control  components. 

V  =  Instantaneous  linear  velocity  along  x  axis 

W  =  Instantaneous  linear  velocity  along  z  axis 

Q  =  Instantaneous  angular  velocity  along  y  axis 

A  =  instantaneous  angle  of  attack 

®  =  Instantaneous  pitch  angle 
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r  =   Instantaneous  glide  angle 

6  =  Elevator  deflection  perturbation 

e  ^ 

6  =  Thrust  perturbation 

2.   Assumptions  Used  in  Obtaining  Equations  of  Motion 

a.  The  airframe  is  a  rigid  body  (no  aeroelastic 
deflection  of  the  airframe) . 

b.  The  earth  is  fixed  in  space,  and  the  earth's 
atmosphere  is  fixed  with  respect  to  the  earth. 

c.  The  mass  of  the  airplane  remains  constant 
during  any  particular  dynamic  analysis. 

d.  The  x-z  plane  (vertical  plane)  is  a  plane  of 
symmetry. 

e.  Disturbances  from  steady  flight  are  small 
enough  so  that  the  products  and  squares  of  changes  in  ve- 
locities are  negligible  in  comparison  to  the  changes  them- 
selves, and  the  disturbance  angles  are  small  enough  so  that 
the  sines  are  equal  to  the  angles  in  radians  and  the  cosines 
are  equal  to  one.   Products  of  these  angles  are  also  approx- 
imately zero.   Because  these  disturbances  are  small,  the 
change  in  air  density  encountered  by  the  airplane  during 
any  disturbance  is  considered  to  be  zero. 

f.  During  steady  flight  conditions,  the  airplane 
is  assumed  to  be  flying  with  wings  level  and  all  components 
of  velocity  zero  except  for  the  linear  velocity  along  the  x 
axis.   (Equilibrium  flight) 
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g.   The  flow  is  quasi-steady.   (The  air  flow  pat- 
tern around  the  airframe  instantaneously  changes  in  a 
steady  flow  pattern  as  the  airplane  changes  its  orienta- 
tion with  respect  to  its  flight  path.) 

3.   The  Longitudinal  Equations  of  Motion 


v(t)  =  Xvv(t)  +  X  q(t)  -  g6(t)  +  X^OJft) 


(3.1) 


+  X-to(t)  +  X.  6  (t)  +  X6  6  (t) 
e  T 


03(t)  =  Zvv(t)  +  VQq(t)  +  Z  q(t)  +  Z^O^t) 


+   z^w(t)  +  z6  6e(t)  +  z6  6  (t) 
e  T 


q(t)  =  Mvv(t)  +  M  q(t)  +  J^OJ(t)  +  Mj,05(t) 


+  M6  8  (t)  +  m6  6  (t) 

e  T 


where  the  subscripted  X,  Z,  and  M  are  the  Dimensional  Sta- 
bility Derivatives  of  the  airplane  which  are  parameters  for 
a  specific  airplane  in  a  specific  flight  regime  (e.g., 
landing  regime,  subsonic  regime,  etc.). 
By  substituting  the  relationships 


(j)(t)    =   v  a(t) 
q(t)  =  6(t) 


°  (3.2) 
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the  equations  of  motion  can  be  rewritten  in  the  more  fa- 
miliar terms  of   v,  OL,    and   0   which  follow: 


v(t)  =  xvv(t)  -  g6(t)  +  X  6(t)  +  xaa(t) 


+  x-a(t)  +  x.  6  (t)  +  x,  6m(t) 
e  T 


z  z        z 

d(t)  =  —  v(t)  +  0(t)  +  —  8(t)  +  —  a(t) 
o                o         o 

Zd  (1  Z6e       Z*T                                                  (3'3) 

+  —  &(t)  +  — -  6   +  — i  6m 

v  v  e    v    T 

o  o        o 


9(t)  =  Mvv(t)  +  M  0(t)  +  MaOi(t)     +   M&Q!(t) 


+  M6  6  (t)  +  Mg  6  (t) 
e  T 


In  a  rigorous  mathematical  sense,  these  longitudinal 
equations  of  motion,  which  deal  with  perturbations  from 
equilibrium  flight,  are  applicable  only  to  infinitesimal 
disturbances;   however,  aerodynamic  experience  has  shown 
that  quite  accurate  results  can  be  obtained  by  applying 
these  equations  to  disturbances  of  finite,  non-zero 
magnitude. 

Since  this  problem  deals  with  the  landing  approach, 
two  additional  parameters  must  be  considered.   One,  the 
altitude  of  the  airplane,  is  of  primary  importance,  and  an 
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equation  relating  actual  airplane  altitude  to  the  flight 
dynamics  is  required.   Such  a  relationship  can  be  approx- 
imated by 

h(t)  =  vQr(t)  (3.4) 

where  h  is  defined  as  the  instantaneous  vertical  distance 
above  the  ground  of  the  aircraft  wheels.   By  definition 
(see  Figure  4) 

T(t)  =  y     +  y(t) 

°  (3.5) 

y(t)  =  6(t)  -  a(t) 

and  therefore, 

h(t)   =  vQyo  +  vo(6(t)  -  a(t))  (3.6) 

The  second  parameter  is  ground  effect  —  defined  as  the 
effect  on  the  airplane  dynamics  as  the  airplane  nears  the 
vicinity  of  the  stationary  plane  (the  ground).   This  ef- 
fect on  an  airplane  is  non-linear  and  is  dependent  on  the 
size  and  shape  of  the  specific  airplane. 

The  formulation  of  a  general  mathematical  model  to 
describe  the  motions  of  any  airplane  in  the  landing  ap- 
proach is  now  complete.   All  that  remains  is  application 
to  a  specific  airplane. 

B.    THE  SPECIFIC  MATHEMATICAL  MODEL 

The  main  problem  associated  with  the  use  of  the  gener- 
al mathematical  model  is  determination  of  the  numerical 
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values  for  the  Dimensional  Stability  Derivatives  for  the 
specific  airplane.   Because  of  the  vast  amount  of  experi- 
mental data  available  on  the  airplane,  both  from  wind  tun- 
nel studies  and  actual  flight  testing,  the  McDonald 
Douglas  F-4J  Phantom  II  jet  fighter  was  selected.   Appendix 
A  contains  representative  data  for  this  airplane  in  the 
landing  configuration. 

As  the  F-4J  Dimensional  Stability  Derivatives,  X  ,  X- , 

q   Of 

Z  ,  Z*,  and  M*   are  zero  and  ground  effects  are  neglected, 
the  longitudinal  equations  of  motion  and  the  height  equa- 
tion can  be  simplified  as  follows: 


v(t)  =  Xvv(t)  +  Xa«(t)  -  g6(t)  +  X6  6T(t) 

a(t)  =  ~  v(t)  +  ^2  a(t)  +  6(t)  +  -^  6e(t)  +  v^T(t) 
o         o  o  o 


(3.7) 


6(t)  =  M  v(t)  +  M  <*(t)  +  M-a(t)  +  M  0(t)  +  M,  6  (t) 
v         0!         Of         3         6   e 


e 


h(t)  =  v  (8(t)  -  a(t) )  +  v  y 

\     i  Q  \   \  /      \  /  /      o  o 


Defining 


x  =    ©  (3.8) 
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and 

/  6e\ 

(3.9) 

the  equations  are  manipulated  into  the  state  variable  form 
x(t)  =  Ax(t)  +  Bu(t)  +  C  (3.10) 

with  the  result 


x(t)  =  x(t) 


all 

a12 

al3 

0 

0 

a21 

a21 

0 

1 

0 

0 

0 

0 

1 

0 

a41 

a42 

0 

a44 

0 

v° 

a52 

a53 

0 

0 

o   b12 

b21    b22 
+   I  |  u(t)   + 


b41    b42 


0 


(3.11) 
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where 


a,,  =  X 
11    v 


52       o 


a12  "  xa 


l53 


v 


o 


a13  =  "5 


v 


21    v 


b12  "  X6r 

ze' 


21    v 


a22  "  v 
o 

M-Z 
&   v 

a.,  =  M  +  

41     v    v 


a42  ■  M<*  +  -^T- 


22    v 


b4i  ■  Mfi  +  -r-5 

e     o 


MaZ6 


42     v 


a.-  =  M   +  Mi. 
43     q     tt 


ccl  =  v  y 
51    o  o 


Appendix  B  contains  the  numerical  values  for  these  elements. 
Heretofore  in  the  formulation,  it  was  assumed  that  the 
two  perturbation  controls,  elevator  deflection  and  thrust, 
are  instantaneously  available  to  control  the  airplane.   In 
reality,  both  controls  have  associated  dynamics,  which  are 
aerodynamically  referred  to  as  elevator  actuator  lag  and 
lag  between  the  thrust  command  and  thrust.   These  effects 
can  be  approximated  by  use  of  the  following  differential 
equations : 


6e(t)  =  -f-    (6   (t)  -  6  (t)) 

6     c 
e 

6T(t)  =  Y*~    <6T  (t)  "  6T(t)) 
6T 


(3.12) 
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where  the  subscripted  T's  are  the  appropriate  component 

time  constants  and  6   and  6_   are  the  elevator  deflection 

e       T 

c        c 

and  thrust  commands  respectively.  In  this  investigation 
it  will  be  assumed  that  the  controls  are  instantaneously 
available;  i.e.,  the  dynamics  will  be  neglected. 

This  completes  the  mathematical  modeling  of  the  F-4J 
airplane  in  the  landing  configuration.   The  model  was  test- 
ed on  an  IBM  System/360  digital  computer,  using  a  Fourth- 
Order  Runge-Kutta  routine  to  solve  the  differential  equa- 
tions.  The  responses  resulting  from  various  control  in- 
puts were  consistent  with  actual  airplane  responses. 
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IV.   OPTIMAL  CONTROL  THEORY 

This  chapter  will  not  attempt  to  present  the  entire 
field  of  optimal  control  theory,  but  rather  will  assume 
that  the  reader  is  familiar  with  it;   hence,  only  those 
techniques  pertinent  to  this  investigation  will  be  dis- 
cussed.  Some  definitions  will  first  be  presented,  follow- 
ed by  mathematical  techniques  employed  to  obtain  the  op- 
timal control. 

A.    OPTIMAL  CONTROL  DEFINITIONS 

1.  Control  History 

The  history  of  control  input  values  during  the 
time  interval  [t  , tf ]  is  denoted  by  u. 

2.  State   Trajectory 

The  history  of  state  values  during  the  time  in- 
terval [t  , t _]  is  denoted  by  x. 
of  — 

3 .  Admissible  Control 

U  denotes  the  set  of  all  control  histories  which 

satisfy  the  physical  control  constraints  during  the  time 

interval  [t  ,  t ,_]  . 
o   f 

4.  Admissible  Trajectory 

X  denotes  the  set  of  all  state  trajectories  which 

satisfy  the  physical  state  variable  constraints  during  the 

time  interval  [t  ,  t ,.]  . 

o   f 
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5.  The  Performance  Measure 

J  denotes  a  scalar  measure  of  the  performance  of 
a  system  when  a  control  history  is  applied. 

6.  Optimal  Control 

A  control  which  causes  the  system  to  follow  an 
admissible  trajectory  which  minimizes  the  performance  mea- 
sure is  called  an  optimal  control  and  is  denoted  by  u*. 

7 .  Optimal  Trajectory 

The  admissible  trajectory  which  results  when  an 
optimal  control  is  applied  to  the  system  is  called  an  op- 
timal trajectory  and  is  denoted  by  x*. 

8.  Tracking  Problem 

A  problem  wherein  the  intent  is  to  maintain  the 
state  trajectory  as  close  as  possible  to  a  desired  trajec- 
tory —  denoted  by  r  —  in  the  interval  [t  , t-]  is  called 
a  tracking  problem. 

B.    THE  LINEAR  TRACKING  PROBLEM 

Since  the  automatic  landing  problem  as  formulated  in 
Chapter  II  is  directed  toward  controlling  the  airplane  in 
a  desired  manner  during  the  final  phases  of  a  landing,  and 
the  mathematical  model  of  the  airplane  as  presented  in 
Chapter  III  is  linear,  the  problem  can  be  considered  as  a 
linear  tracking  problem. 
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It  has  been  shown  (Ref .  6  and  7) ,  where  complete  deri- 
vations can  be  found,  that,  given  a  set  of  state  equations 


x(t)  =  A(t)  x(t)  +  B(t)u(t) 


(4.1) 


where 

x(t)  is  the  state  vector  of  dimension  n 
u(t)  is  the  control  vector  of  dimension  m 
A(t)  is  a   nxn  matrix 
B(t)  is  a   nxm  matrix 
and  the  performance  measure 


J  =  \   ||x(tf)  -  r(tf)  || 

H 


where 


llx(T)  -  r(T) 


2  2 

x(t)  -  r(t)  ||      +  ||u(t)|| 

Q(t)  R(t) 


dt 


(4.2) 


'Q(T) 


nT 


x(T)  -  r(T) 


Q(T) 


x(T) 


-  r(T)]  , 


the  final  time  t^-  is  fixed,  ^(^f)  is  free,  u(t)  is  uncon- 
strained, H  and  Q(t)  are  real  symmetric  positive  semi- 
definite  matrices,  and  R(t)  is  a  real  symmetric  positive 
definite  matrix,  that  the  optimal  control  exists  and  is 
unique.   The  optimal  control  is  given  by 


u*(t)  =  F(t)  x*(t)  +  g(t) 


(4.3) 
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where  F(t)  is  the   mxn  matrix  of  feedback  gains,  and  g(t) 
is  the   mxl   command  signal  vector  which  is  dependent  on 
the  system  parameters  and  the  desired  state  trajectory. 
Figure  5  is  a  block  diagram  of  the  plant  and  optimal  con- 
troller. 
F(t)  and  g(t)  are  given  by 

F(t)  =  -  R(t)"1  BT(t)K(t)  (4.4) 

g(t)  =  -  R(t)"1  BT(t)  s(t)  (4.5) 

where  K(t)  is  the  solution  of  the  Riccati-type  matrix  dif- 
ferential equation 


K(t)  =  -  K(t)A(t)  -  AT(t)K(t)  -  Q(t) 


+  K(t)B(t)R  1(t)BT(t)K(t) 


(4.6) 


with  boundary  conditions 

K(tf)  =  H  (4.7) 

and  s_(t)  is  the  solution  of  the  linear  vector  differential 
equation 

s(t)  =  -  [AT(t)-K(t)B(t)R_1(t)BT(t)]s(t)+Q(t)r(t) 

(4.8) 
with  boundary  conditions 

s(tf)  =  -  Hr(tf)  (4.9) 
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PLANT 


U 


CONTROLLER 


jk 


L 


F(t) 


FIGURE   5 
BLOCK  DIAGRAM    OF  PLANT  AND 

CONTROLLER 
IN    A   LINEAR  TRACKING  PROBLEM 
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Since  the  boundary  conditions  at  the  final  time  tf  are 
known,  these  two  matrix  differential  equations  must  be  in- 
tegrated backwards  in  time  from  t   to  t  . 

f     o 

It  was  assumed  that  all  states  of  the  system  can  be 
measured. 


C.    MODIFICATIONS  TO  SPECIFY  STATE  PLANT 

The  previous  discussion  indicated  that  the  optimal  con- 
trol in  a  linear  tracking  problem  exists  and  can  be  obtain- 
ed, assuming  that  the  plant  is  of  the  form 

x(t)  =  A(t)x(t)  +  B(t)u(t)  (4.10) 

Since  the  plant  model,  as  presented  in  Chapter  III,  is  in 
the  form 

x(t)  +  Ax(t)  +  Bu(t)  +  c  (4.11) 

some  modifications  are  necessary  in  order  to  adapt  the 
specific  problem  to  the  theory  presented. 

This  modification  can  be  accomplished  by  observing 
that  the  c_  vector  contains  only  one  non-zero  element,  i.e., 
the  differential  equation  describing  the  altitude  of  the 
airplane 

h(t)  =  vQyo  +  vo(e(t)  -  a(t))  (4.12) 

contains  the  constant  term  v  y 

o  o 
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By  defining 

h(t)  =  h(t)  -  h  (t) 


(4.13) 


where 


h  (t)  =  h(0)  +  v  y   t 
e  o  o 


(4.14) 


is  the  nominal  equilibrium  altitude  of  the  airplane  during 


the  interval  Ct  ,  t ,.]  .   It  follows  that 


h  (t)  =  v  y 
ev  '     o  o 


h(t)   =  h(t)  -  he(t) 


h(t)   =  v  y     +   v  (0(t)  -  a(t))  -  v  y 
v  '  o  o    o  o  o 


h(t)   =  vo(9(t)  -  a(t)) 


(4.15) 
(4.16) 
(4.17) 
(4.18) 


So  h(t),  as  defined,  is  the  perturbation  of  the  altitude  of 
the  airplane  about  the  nominal  equilibrium  altitude. 

With  this  modification,  a  new  set  of  states,  hereafter 
referred  to  as  the  revised  states  (to  distinguish  them  from 
the  actual  states) ,  can  be  defined  as 


x   = 


x 


(4.19 
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and  the  revised  state  equations  become 


x(t)  =  Ax(t)  +  Bu(t) 


(4.20) 


where  the  elements  of  A  and  B  remain  identical  to  those  pre- 
viously formulated  in  Chapter  III,  equation  (3.11).   This 
revised  plant  is  now  in  the  proper  form  for  application  of 
the  optimal  control  theory.   The  solution  then  becomes 


u*(t)  =  P(t)x*(t)  +  g(t) 


(4.21) 


from  which  the  behavior  of  the  actual  states  can  be  deter- 
mined by  substituting  the  definition  of  x(t) .   This  substi- 
tution produces  the  optimal  control 


u*(t)  =  F(t)x*(t)  -  F(t) 


+  g(t) 


(4.22) 


By  this  simple  mathematical  manipulation,  the  problem  can 
be  solved  by  applying  optimal  control  theory  to  the  revised 
state  equations  and  then  expressing  the  results  in  terms  of 
the  actual  states. 
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V.   SPECIFICATIONS,  DESIRED  TRAJECTORIES, 
PERFORMANCE  MEASURE 

In  this  chapter  the  applicable  specifications  are  in- 
troduced and  a  set  of  desired  state  trajectories  for  the 
problem  are  derived.   Many  of  the  results  presented  are 
based  on  the  author's  experience  as  a  pilot  operating  the 
F-4  and  other  jet  airplanes.   The  final  section  of  this 
chapter  discusses  the  significance  of  the  performance  meas- 
ure in  obtaining  realistic  results. 

A.    PROBLEM  SPECIFICATIONS 

A  successful  automatic  landing  requires  that  certain 
conditions  relative  to  the  airplane  and  its  environment  be 
satisfied.   These  conditions,  which  assume  the  role  of  spe- 
cifications for  the  problem,  are  formulated  in  terms  of  per- 
formance requirements  and  constraints  on  the  system  states 
and  controls,  together  with  specifications  regarding  related 
parameters.   The  following  specifications  are  of  primary 
importance  to  the  problem  under  consideration. 

1.   Velocity 

The  instantaneous  velocity  of  the  airplane  during 
the  landing  portion  of  flight  must  remain  above  1.1  vstaii 
or  V(t)  ^  214.5  ft/sec.   The  upper  velocity  limit  depends 
on  the  structural  limit  with  landing  gear  and  flaps  extend- 
ed;  flight  experience  indicates  that  a  reasonable  upper 
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limit  is  1.1  v   or  V(t)  ^  245  ft/sec.   In  terms  of  the 
specified  state,  v(t),  which  is  the  perturbation  velocity 
about  equilibrium  flight,  these  limitations  are 

-8.5  ft/sec  ^  v(t)  £  22  ft/sec   ,   tc[t  ,  tf  ]  . 

2.  Angle  of  Attack 

The  angle  of  attack  must  remain  below  0.9  A    ,, 

stall 

or  A(t)  =s  0.44  radians.   No  stringent  lower  limit  exists 
but,  as  in  the  velocity  case,  experience  indicates  that  a 
reasonable  minimum  value  is  approximately  0.750*  or 
A.(t)  ^0.25  radians.   As  a  result,  the  perturbation  angle 
of  attack  limits  were  established  as 

I  a(t)  |   <:  0.11  rad   ,    tcCt  ,  t-]  . 

3 .  Pitch  Angle 

During  the  landing  portion  of  flight,  the  instan- 
taneous pitch  angle  is  closely  related  to  the  angle  of 
attack  but  does  not  have  stringent  limitations  on  its  ex- 
cursions.  However,  at  the  actual  touchdown  point,  limita- 
tions do  exist  primarily  to  prevent  the  airplane  from 
either  landing  nose  wheel  first  (low  pitch  angle)  or  tail 
first  (high  pitch  angle) .   For  the  F-4J  airplane  these  lim- 
itations, in  terms  of  the  perturbation  pitch  angle  from 
equilibrium  flight  conditions,  were  established  as 

-0.20  rad  s  9(t)  <  0.25  rad   ,    te[t  ,  tf]  . 
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4.   Pitch-Angle  Rate 

Flight  experience  has  indicated  that  this  para- 
meter should  be  held  to  a  minimum  for  pilot  comfort.   (See 
discussion  under  Related  Specifications,  below.)   Realistic 
limits  for  the  pitch  angle  rate  were  established  as 

I  © ( t)  |  £  0.08  rad/sec   ,    te[t  ,t  ]  . 

3.   Altitude 

In  the  problem  formulation,  Chapter  II,  it  was 
assumed  that  the  altitude  of  the  airplane  at  the  initial 
time  was  within  the  prescribed  Category  II  ILS  "window"  or 

88  ft  ^  h(t  )  <  112  ft 
o 

In  terms  of  perturbations  about  equilibrium  flight  condi- 
tions, this  limitation  can  be  expressed  as 

|  h(t  )|  ^  12  ft 
v  o 

As  the  airplane  approaches  the  actual  touchdown  point,  how- 
ever, altitude  excursions  about  desired  conditions  should 
satisfy  more  stringent  limits.   As  a  result,  the  following 
limitations  were  established  for  altitude  perturbations 
about  a  desired  altitude  trajectory,  in  the  intervals  in- 
dicated.  (See  Section  B,  this  chapter,  for  desired  trajec- 
tory. ) 

In  the  interval  from  the  decision  altitude  to  the 

flare  initiation  altitude 

|h(t)  -  hd(t)|  <  12ft  ,   t€[tQf  t^], 
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and  in  the  interval  from  the  flare  initiation  altitude  to 
actual  touchdown 

I  h(t)  -  hd(t)|  £  5ft,     t€[t1#tf]  , 

where  h,  (t)  is  the  desired  revised  altitude  state  trajec- 
tory, and  t,  is  the  time  at  which  the  flare  phase  begins. 

6 .  Elevator  Deflection 

Physical  limits  for  elevator  travel  exist,  and  for 
the  F-4J  airplane  in  equilibrium  flight  in  a  landing  ap- 
proach, were  established  as 

-0.26  rad  ^  &e(t)  s  0.22  rad,     tc[tQ,tf]. 

7.  Thrust 

Physical  limits  on  thrust  available  exist  for  the 
F-4J  airplane  and  are  dependent  on  the  selection  of  power; 
i.e.,  MILITARY  engine  operation  or  AFTER  BURNER  (A/B)  en- 
gine operation.   Flight  experience  has  shown  that  realistic 
limits  for  thrust  perturbations  about  equilibrium  flight 
conditions  in  the  landing  portion  of  flight  can  be  estab- 
lished as 

I  5  (t)|  <  3000  lbs,    te[t  ,tf]  . 

This  range  is  well  within  the  thrust-available  range  of 
the  F-4J  airplane  without  A/B  selection. 
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8.   Related  Specifications 

Several  related  performance  requirements  and  con- 
straints are  applicable  to  the  problem  and  are  presented 
in  the  following  listing: 

a.  The  actual  touchdown  point  must  be  within  150 
feet  of  the  desired  touchdown  point  on  the  landing  surface. 
Since  the  problem  was  formulated  in  the  time  domain,  this 
restriction  was  established  by  requiring  that  the  actual 
touchdown  occur  within  ±0.65  seconds  of  the  desired  touch- 
down time . 

b.  At  the  actual  touchdown  time  the  airplane  must 
have  no  tendency  to  float,  or,  in  other  words,  the  airplane 
should  have  a  positive  rate  of  descent.   In  addition,  the 
rate  of  descent  at  touchdown  must  be  within  the  structural 
sink  rate  limits  for  the  airplane  (Appendix  A) .   Experience 
indicates  that  realistic  rate  of  descent  limits  for  the  F-4J 
at  the  actual  touchdown  time  are 

-9ft/sec  <  h(tf)  <  -3ft/sec 

c.  The  normal  acceleration   of  the  airplane  -- 
defined  as  a  perturbation  from  the  gravitational  accelera- 
tion of  equilibrium  flight  (l.Og)  --  can  be  approximated  by 

n(t)  =  -^(e(t)  -  ■(*)) 


l 
The  airplane  acceleration  perpendicular  to  the  hori- 
zontal reference  plane. 
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and  has  physical  limits  based  on  the  structural  strength  of 
the  airplane.   These  limits  are  normally  not  attained  in 
the  landing  portion  of  flight,  due  to  the  low  velocity  and 
control  effectiveness  in  this  regime.   However,  the  func- 
tion n(t)  gives  an  indication  of  the  smoothness  of  the  land- 
ing, and  its  effect  can  actually  be  felt  by  the  pilot  during 
the  landing:   therefore,  to  maintain  the  normal  accelera- 
tion as  closely  as  possible  to  the  equilibrium  flight  condi- 
tion of  l.Og  the  following  limitation  was  established 

I  n(t)|   £  0.2g,     t([t,tf]  . 

B.    DESIRED  TRAJECTORY 

The  problem  has  been  formulated  as  a  tracking  problem; 
hence  the  desired  state  trajectory  must  be  determined.  This 
section  presents  a  general  desired  trajectory  appropriate 
for  any  airplane  in  the  landing  portion  of  flight,  and  in- 
dicates how  these  trajectories  are  applied  to  the  F-4J 
landing  problem. 

1.   General  Desired  Trajectory 

Since  the  altitude  of  the  airplane  is  of  great 
importance  during  the  landing,  it  is  an  obvious  starting 
point  in  the  formulation  of  a  set  of  desired  state  trajec- 
tories.  In  Chapter  II  it  was  stated  that  this  investiga- 
tion would  consider  only  that  portion  of  the  landing  from 
the  decision  altitude  (100  feet)  to  the  touchdown  point. 
Referring  to  Figure  2,  this  portion  of  the  landing  includes 
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a  section  of  the  capture-track  phase,  the  flare  phase,  and 
the  land  phase.   In  Chapter  III  it  was  shown  that  the  pro- 
blem could  be  formulated  in  the  time  domain.   These  two 
considerations  lead  to  a  general  altitude  trajectory  (spe- 
cified in  the  time  domain)  consisting  of  a  constant  descent 

from  the  decision  altitude  at  t  =  t   =  0  to  a  selected 

o 

flare  point  at  t  =  t, ,  followed  by  a  flare  to  touchdown  at 
t  =  tf. 

In  the  optimal  control  theory  presented  in  Chapter  IV, 
one  of  the  necessary  requirements  was  that  the  final  time 
(tr.)    be  fixed.   To  alleviate  the  obvious  disadvantage  that 
this  requirement  imposes,  it  seemed  appropriate  to  extend 
the  final  time  beyond  the  actual  touchdown  time  to  ensure 
that  a  landing  occurs  prior  to  t^.   By  so  doing,  an  imag- 
inary flare  plane  below  the  actual  landing  plane  was 
established.   This  implies  that  the  airplane  nominally  will 
land  before  the  final  time  specified.   Figure  6  depicts 
this  general  concept  in  graphical  form,  where  the  capture- 
track  phase  is  the  interval  [0,t, ),  the  actual  flare  phase 
is  the  interval  [t,  ,t9),  and  the  actual  land  phase  is  t„. 
The  imaginary  flare  phase  is  the  interval  [t,,t  )  and  the 
imaginary  land  phase  is  at  the  time  t... 

The  rate  of  altitude  change  (rate  of  descent),  h(t), 
must  be  considered  in  conjunction  with  the  general  altitude 
trajectory,  even  though  it  is  not  a  state  variable.   Ideal- 
ly, during  the  capture-track  phase  [0,t,  ),  the  rate  of 
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descent  will  be  a  constant,  while  during  the  flare  phase, 
it  will  steadily  decrease  so  that  at  the  actual  touchdown 
time  (t2),  the  rate  of  descent  will  be  within  the  estab- 
lished limits.   Utilizing  this  consideration,  the  set  of 
equations 


hd(t)  =  he(t)  +  h(0)  +  vQyot  (5.1) 

hd(t)  =  he(t)  =  vQyo  (5.2) 

which  describe  the  desired  altitude  trajectory,  h,(t),  and 
rate  of  descent,  h,  (t),  for  the  interval  l_0,t,)  were  for- 
mulated. 

For  the  interval  [t, ,t  ]  ,  the  desired  altitude  and 
rate  of  descent  were  defined  as 

hd(t)  =  h(0)  +  voyot  +  -1  (e  1  1   -l)  -L2(t-tl) 

(5.3) 


Since  the  revised  altitude  state  was  defined  in  Chapter  IV 
to  be 

h(t)  =  h(t)  -  he(t)  (5.5) 

equations  (5.1)  and  (5.2)  became 

hd(t)  =  hd(t)  =0,     t€[0,t),  (5.6) 
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Figure  6 
General  Desired  Altitude  Trajectory 


Plane 


and 


L9  •  L,  ( t-t, )  ^ 

hd(t)  =  iKe  -v-  ^ct-^j 


(5.7) 


h,(t)  =  L. 


L,  (t-t,) 


^e      x  -lj  ,    teCfc^tf] 


(5.8) 


Equations  (5.6)  and  (5.7)  represent  the  desired  revised 
altitude  trajectory  in  the  interval  [0,tf]  in  terms  of  per- 
turbations about  the  equilibrium  altitude. 
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These  general  expressions  for  the  desired  revised  rate 
of  descent,  h,(t),  provide  insight  into  the  form  of  the 
desired  angle  of  attack  trajectory,  a,,  and  pitch  angle 
trajectory,  &    ,  since  from  equation  (4.18) 

h(t)  =  vo(6(t)  -  tt(t))  .  (5.9) 

Therefore, 

ed(t)  -  ad(t)  =  0,    tcCO,^)  (5.10) 


L2  f   Li(t"ti) 


and 


6d(t)  -  ad(t)  =  rf-  {e   x      -1/.   tcCt1,tf].  (5.11) 

o 

By  definition  (Chapter  II)  ,  Oi   and  8   are  zero  when  the  air- 
plane is  in  equilibrium  flight,  hence  equation  (5.10)  is 
consistent.   The  problem  remains  to  find  equations  for  Oi 
and  0  in  the  interval  [t,  ,  t^.]  which  are  consistent  with 
the  requirements  of  equation  (5.11).   Before  proceeding, 
recall  that  the  desired  pitch-angle-rate  trajectory,  8 
must  also  be  considered.   An  acceptable  expression  for  the 
desired  pitch-rate  trajectory  is  a  linear  function  during 
the  flare  phase  interval  [t,,tf]  given  by 

6d(t)  =  L3(t-tL)  (5.12) 
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Integrating  equation  (5.12)  and  applying  the  desired  bound- 
ary condition  at   t  =  t, ,  i.e.,   0  (t.)  =  0 

9d(t)  =  ~T    (t  "  tl)2#  (5.13) 

Substituting  equation  (5.13)  into  (5.11)  yields 

L  2    L    L  (t-t  ) 

«d(t)  =  -J-   (t  -  t  P  -  -f-  [e  -l)  (5.14) 

o 

as  the  expression  for  the  desired  angle  of  attack  in  the 
interval  [t,,tf]. 

The  final  state  component  trajectory  which  must  be 
considered  is  the  velocity  trajectory.   Experience  from 
both  an  operational  and  a  safety  standpoint  has  establish- 
ed that  the  airspeed  during  a  landing  should  remain  es- 
sentially constant.   This  implies  that  v.,  the  desired 
velocity  trajectory,  should  be  zero  during  the  entire  in- 
terval of  interest,  hence 


vd(t)  =  0,    te[0,tf]  .  (5.15) 


2 .   Specific  Desired  Trajectories 

In  the  preceding  section  a  general  set  of  desired 
state  trajectories  —  applicable  to  any  airplane  —  were 
formulated.   What  remains  is  application  to  the  F-4J  air- 
plane so  that  the  unknown  constant  terms  can  be  evaluated. 
Since,  in  the  interval  [o,t,),  the  desired  state  trajec- 
tory is  zero,  the  following  discussion  will  be  concerned 
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with  the  interval  [t,,t  ]  where  the  equations  of  interest 
are  (5.3),  (5.4),  (5.7),  (5.8),  (5.12),  (5.13)  and  (5.14) 
—  hereafter  referred  to  as  the  reference  equations. 

The  first  consideration,  however,  was  to  determine  a 
realistic  time  frame  for  the  problem.   The  equilibrium  air- 
speed, v  ,  of  the  F-4J,  as  given  in  Appendix  A,  is  223 
ft/sec.   Assuming  that  the  airplane  is  at  the  nominal  de- 
cision altitude  (100  feet)  at  t  =  0  and  that  the  airplane 
remained  in  equilibrium  flight  on  a  negative  three-degree 
glide  slope,  the  actual  time  of  touchdown  from  equation 
(4.14)  is  approximately  8.57  seconds.   This  assumes  that 
there  is  no  flare  phase.   Although  the  standard  operating 
procedure  in  the  U.S.  Navy  is  to  operate  the  F-4J  airplane 
with  no  flare,  this  investigation  will  include  a  flare 
which  commences  at  t,  =6  seconds;   at  this  time  the  air- 
plane is  at  a  nominal  altitude  of  approximately  30  feet. 
It  follows  then  that  a  reasonable  time  frame  for  the  prob- 
lem is  10  seconds,  and  that  the  interval  [t,,tf]  is  4  sec- 
onds in  duration. 

Next,  values  for  the  unknown  constants  in  the  refer- 
ence equations  must  be  obtained.   Boundary  conditions  for 
the  reference  equations  at  t  =  t,  =  6  seconds  are  known, 
and  final- time  boundary  conditions  can  be  selected  to  pro- 
vide realistic  trajectories.   However,  recalling  that  the 
actual  touchdown  of  the  airplane  occurs  prior  to  the  final 
time,  and  that  the  specifications  refer  to  this  actual 
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touchdown  time,  another  set  of  boundary  conditions  at 
t  =  tp  must  be  satisfied.   The  conditions  at  t  =  L  are  not 
exactly  given,  but  realistic  values  can  be  selected  to  meet 
the  necessary  specifications.   As  a  result,  values  for  the 
unknown  constants  L,  ,  L2,  and  L-.  of  the  reference  equations 
were  found  by  a  trial-and-error  method.   Values  of  these 
constants  were  arbitrarily  chosen  to  satisfy  boundary  con- 
ditions at  t  =  t,  and  then  substituted  into  the  reference 
equations  to  obtain  the  actual  touchdown  time,  t~,  and  the 
corresponding  values  of  the  parameters  at  t  =  L.   This 
process  was  continued  until  values  of  L, ,  L?,  and  L_  were 
found,  which  resulted  in  realistic  trajectories.   The  values 

L,  =  0.25 
L2  =  4.68 
L3  =  0.016 

which  were  obtained  result  in  a  desired  actual  touchdown 
time,  (t?),  of  9.3  seconds.   The  reference  equations  which 
apply  to  the  interval  [6,10]  are  then 

hd(t)  =  h(0)  +  vQyot  +  £|§  (e*25(t-6)-l)  -4.68(t-6) 

(5.16) 
hd(t)  =  vQyo  +  4.68  (e' 25 ( t_6) -l)  (5.17) 
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Hd(t)  =  ftff  (e-25(t-6)-l)-4.68(t-6)  (5.18) 

Ed(t)  =  4.68  (e*25(t_6)-l)  (5.19) 


6d(t)  =  0.016(t-6)  (5.20) 

0d(t)  =  0.008(t-6)2  (5.21) 

ad(t)  =  0.008(t-6)2  -  ^^  (V25(t_6)-l)  ,      (5.22) 

o 

while  in  the  interval  [0,6)  the  applicable  equations  are 


hd(t)  =  h(0)  +  VQyot  (5.23) 


hd(t)  =  v  y  (5.24) 


hd(t)  =  hd(t)  =  9d(t)  =  9d(t)  =  ad(t)  =  0       (5.25) 


and 


vd(t)  =  0,    tc[0#10]  .  (5.26) 

Table  A  summarizes  the  solutions  of  equations  (5.16)  through 
(5.22)  at  the  flare  initiation  time  (t, ),  the  actual  touch- 
down time  (t2),  and  the  final  time  (tf). 
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Table  A 

Solutions  to  Reference  Equations  at  Flare  Initiation 
Time,  Actual  Touchdown  Time,  and  Final  Time 

Solutions  at 


Parameter 

Eqi 

uation   No. 

t=6 

sec 

t=9.3    sec 

t=10    sec 

hd(ft) 

5.16 

29.947 

0 

-3.286 

• 

h , (ft/sec) 

5.17 

-11.676 

-5.667 

-3.621 

hd(ft) 

5.18 

0 

8.582 

13.469 

• 

h , (ft/sec) 

5.19 

0 

6.009 

8.055 

• 

0 , (rad/sec) 

5.20 

0 

0.053 

0.064 

©d(rad) 

5.21 

0 

0.087 

0.128 

«d(rad) 

5.22 

0 

0.060 

0.092 

Figure  7  presents  a  plot  of  the  actual  desired  altitude 
trajectory  with  the  equilibrium  altitude  trajectory  includ- 
ed for  comparison.   Figure  8  illustrates  the  revised  desired 
altitude  trajectory  h,,  while  Figure  9  depicts  the  three 
angular  components  of  the  desired  state  trajectory,  a,  , 

0.,,  and  0  .,  and  the  velocity  trajectory  v... 
d       d  d 
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Desired  Altitude  Trajectory 
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Desired  Revised  Altitude  Trajectory 
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Desired  Velocity,  Angle  of  Attack,  Pitch  Angle,  and 
Pitch  Rate  Trajectory 
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C.    PERFORMANCE  MEASURE 

The  mathematical  form  of  the  performance  measure  in- 
troduced in  Chapter  IV  was 

2 

J   =    1/2    ||    x(t    )     -    r(t    )     || 

H 


*i 


f  r  2  2 

[||   x(t)    -   r(t)     ||  +    ||   u(t)  || 

tQ  Q(t)  R(t)J 


dt 


(5.27) 

One  of  the  advantages  of  this  quadratic  performance  measure 
is  that  the  elements  of  the  H,  Q(t)  and  R(t)  matrices, which 
shall  be  called  the  weighting  matrices,  can  be  related  to 
the  design  parameters  of  the  system  and  chosen  to  satisfy 
the  design  objectives.   If  these  weighting  matrices  meet 
the  requirements  established  in  Chapter  IV,  specifically 
H  and  Q(t)  are  positive  semi-definite  and  R(t)  is  positive 
definite,  the  optimal  control  law  can  be  found,  and  is 
unique  for  the  selected  set  of  weighting  matrices. 

Previous  investigations  dealing  with  similar  airplane 
landing  problems  (Refs.  8  and  9)  have  indicated  that  the 
weighting  matrices  must  be  diagonal  matrices  with  positive 
diagonal  entries.   In  other  words,  the  weighting  matrices 
must  have  non-zero  weighting  values  specified  for  each 
state  and  control.   As  a  result,  in  the  initial  trials  it 
was  assumed  that  the  weighting  matrices  are  diagonal  and 
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time  invariant  over  the  interval  of  interest.   Since  the 
form  of  the  weighting  matrices  was  so  selected,  each  dia- 
gonal element  of  the  respective  matrix  could  be  related  to 
a  specific  state  variable  or  control,  and  therefore,  theo- 
retically, could  be  assigned  a  proper  value  to  meet  estab- 
lished specifications.   Unfortunately,  this  selection  of 
weighting  factors  was  not  simply  accomplished.   In  fact, 
it  became  obvious  during  the  investigation  that  the  state 
trajectory  was  very  sensitive  to  changes  in  the  weighting 
matrices.   (See  related  discussion  in  Chapter  VII.) 

The  selection  of  values  for  the  weighting  matrices  was 
carried  out  by  a  trial-and-error  method  wherein  a  set  of 
weighting  matrices  was  selected,  the  optimal  control  law 
computed,  the  optimal  trajectory  calculated,  and  the  re- 
sults compared  with  established  specifications.   In  effect, 
the  final  weighting  matrix  selection  was  based  on  obtain- 
ing a  realistic  optimal  trajectory  which  conforms  to  the 
specifications . 

Appendix  B  contains  the  numerical  values  used  in  the 
weighting  matrices  for  this  investigation. 
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VI.   INVESTIGATION  PROCEDURE 

Since  the  previous  discussion  has  presented  the  pro- 
cedures used  in  modeling  the  plant,  developing  the  desired 
trajectory,  and  related  topics,  only  the  procedure  used  to 
obtain  the  optimal  control  will  be  presented.   The  actual 
cases  investigated  will  also  be  discussed  in  the  last  sec- 
tion of  this  chapter. 

A.    PROCEDURE 

The  following  procedure  was  used  to  obtain  the  op- 
timal control: 

1.  A  set  of  representative  values  were  selected  for 
the  diagonal  elements  of  the  weighting  matrices  (H,  Q,  and 
R)  of  the  performance  measure,  equation  (5.27). 

2.  Equations  (4.6)  and  (4.8)  were  simultaneously  in- 
tegrated from  t  =  t,  to  t  =  0  to  obtain  values  for  K(t) 
and  s_(t)  respectively. 

3.  Using  the  results  of  step  2,  equations  (4.4)  and 
(4.5)  were  solved  to  obtain  values  for  F(t)  and  g(t)  re- 
spectively. 

4.  The  state  equations  (3.11)  were  integrated  from 
t  =  0  to  t  =  t,.,  using  the  results  of  step  3  to  obtain 
u*(t)  of  equation  (4.22). 

5.  The  optimal  control  and  optimal  trajectory  ob- 
tained from  step  4  were  observed  to  ascertain  compliance 
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with  problem  specifications.   If  the  results  were  not  ad- 
missible or  unrealistic,  the  values  of  the  elements  of  the 
weighting  matrices  were  changed  and  steps  2  through  4  were 
repeated. 

To  implement  the  foregoing  procedure,  a  computer 
program  was  written  in  Fortran  IV  for  use  in  the  IBM  Sys- 
tem/360 digital  computer.   System/360  Scientific  Subrou- 
tines were  used  in  the  program  to  accomplish  the  matrix 
algebra.   An  existing  Fourth-Order  Runge-Kutta  subroutine 
developed  at  the  Naval  Postgraduate  School  was  used  for 
the  required  numerical  integration.   The  computer  program, 
including  the  integration  subroutine,  is  presented  after 
the  appendices. 

B.    CASES  INVESTIGATED 

Two  cases  were  investigated  in  the  initial  study. 
Case  I  assumes  that  the  velocity  of  the  airplane  is  con- 
stant and  elevator  deflection  is  the  only  control  avail- 
able;  Case  II  considers  the  plant  and  controls  as  pre- 
viously formulated.   As  a  result,  two  different  state  mo- 
dels were  used. 

1.   Case  I 

Since  it  is  assumed  that  the  velocity  is  constant 
and  that  elevator  deflection  is  the  only  control  available, 
the  state  model,  equation  (3.11)  reduces  to 
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x(t) 


0 
0 
0 

l53 


l44 
0 


x(t) 


and 


(6.1) 


+ 


u(t)       + 


where 


x 


(6.2) 


u        =      6 


(6.3) 
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By  the  same  argument  as  was  presented  in  Chapter  IV,  this 
state  equation  can  be  manipulated  into  the  revised  state 
equation 

•  . 

x(t)  =  Ax(t)  +  Bu(t)  (6.4) 

where 


x  =     .  (6.5) 


and  the  elements  of  A  and  B  remain  identical  to  those  of 
equation  (6.1) . 

2.   Case  II 

The  state  model,  equation  (3.11),  and  the  revised 
state  model,  equation  (4.20),  as  previously  presented,  were 
used  in  Case  II. 

In  order  to  evaluate  the  optimal  trajectory,  three  dif- 
ferent sets  of  initial  conditions  were  selected.   The  first 
assumes  that  the  airplane  is  in  an  ideal  flight  condition, 
or,  in  terms  of  the  revised  state  model,  all  initial  state 
values  are  zero.   The  other  two  sets  were  selected  to  con- 
form to  the  worst  initial  flight  conditions  that  could  be 
reasonably  expected.   The  following  listing  discusses 
these  initial  flight  conditions. 
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a.  HIGH  AND  FAST.   This  set  of  initial  conditions 
assumes  that  the  airplane  was  at  the  upper  altitude  limit 
of  the  ILS  Category  II  "window"  and  that  the  airplane  ve- 
locity is  above  v  ,  the  equilibrium  airspeed.   In  Case  I, 
the  second  assumption  is  not,  strictly  speaking,  applicable 
since  it  is  assumed  that  velocity  was  constant.   However, 
the  condition  was  simulated  by  assuming  that  the  initial 
angle  of  attack,  Oi(t   ),  and  pitch  angle,  0(t  ),  were  below 
their  respective  equilibrium  values.   This  condition  would 
result  in  actual  flight  if  the  airplane  velocity  were 
above  the  equilibrium  flight  velocity.   It  was  further 
assumed  that  the  initial  pitch  angle  rate,  0(t  ),  was  zero. 

b.  LOW  AND  SLOW.   This  set  of  initial  conditions 
assumes  that  the  airplane  is  at  the  lower  altitude  limit  of 
the  ILS  Category  II  "window"  and  that  the  airplane  velocity 
is  below  the  equilibrium  airspeed.   By  the  same  reasoning  as 
before,  the  second  assumption  was  simulated  for  Case  I  by 
assuming  that  the  initial  angle  of  attack  and  pitch  angle 
were  above  their  respective  equilibrium  values.   Again, 

the  initial  pitch  angle  rate  was  assumed  to  be  zero. 

Table  B  presents  a  summary  of  the  cases  investigated 
and  includes  the  initial  values  used  for  each  case. 
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Table  B 
Summary  of  Cases  Investigated 


Initial  Flight 

Case 

Condition 

IA 

IDEAL 

IB 

HIGH  AND  FAST 

IC 

LOW  AND  SLOW 

IIA 

IDEAL 

IIB 

HIGH  AND  FAST 

IIC 

LOW  AND  SLOW 

T 


Applicable      Initial  State 
State  Equations     Conditions 

(6.1,6.4)         (0,0,0,0)T 

(6.1,6.4)  (-0.03,-0.03,0,112.0) 
(6.1,6.4)     (0.03,0.03,0,88.0) 
(3.11,4.20)        (0,0,0,0,0)T 

(3.11,4.20)  (5.0, -0.03, -0.03,0, 112.0 f 

(3.11,4.20)  (-5. 0, 0.03, 0.03,0, 88. 0)T 
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VII.   RESULTS 

A  general  result  is  first  presented,  followed  by  spec- 
ific results,  in  graphical  form,  for  each  of  the  cases  in- 
vestigated.  The  final  section  of  this  chapter  discusses  a 
problem  encountered  during  the  investigation. 

A.    GENERAL  RESULTS 

The  optimal  control  for  an  all-weather  landing  in  the 
F-4J  airplane  was  derived  (see  Specific  Results  below) 
using  the  techniques  presented.   Since  some  simplifying  as- 
sumptions were  made  to  reduce  the  complexity  of  the  prob- 
lem, the  results  are  not  conclusive,  but  serve  to  demon- 
strate that  the  design  of  an  automatic  controller  for  the 
landing  of  an  airplane  is  feasible  by  formulating  the  prob- 
lem as  a  tracking  problem  and  applying  optimal  control 
theory.   By  systematically  eliminating  some  of  the  assump- 
tions made  in  this  study,  a  closer  approximation  to  the 
actual  airplane  landing  problem  can  be  accomplished. 

Actual  implementation  of  an  automatic  controller  to 
provide  the  optimal  control  derived  by  these  techniques 
would  require  that  a  digital  computer  be  placed  aboard  the 
airplane  which  may  impose  an  undersirable  penalty.   If  so, 
the  optimal  control  still  provides  the  goal  for  the  de- 
sign of  any  sub-optimal  controller  and  hence  is  of  great 
value  to  the  control  engineer. 
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B.    SPECIFIC  RESULTS 

Specific  results  for  each  of  the  cases  investigated 
are  presented  in  graphical  form  and  are  summarized  below. 

1.  The  solution  to  equation  (4.6),  the  K(t)  values, 
are  presented  for: 

a.  Case  I  in  Figure  10 

b.  Case  II  in  Figure  21 

2.  The  solution  to  equation  (4.8),  the  s^(t)  values, 
are  presented  for: 

a.  Case  I  in  Figure  11 

b.  Case  II  in  Figure  22 

3.  The  optimal  control  is  depicted  for: 

a.  Case  IA  in  Figure  12 

b.  Case  IB  in  Figure  15 

c.  Case  IC  in  Figure  18 

d.  Case  IIA  in  Figure  23 

e.  Case  IIB  in  Figure  26 

f.  Case  IIC  in  Figure  29 

4.  The  optimal  and  desired  altitude  trajectories  are 
presented  for: 

a.  Case  IA  in  Figure  13 

b.  Case  IB  in  Figure  16 

c.  Case  IC  in  Figure  19 

d.  Case  IIA  in  Figure  24 

e.  Case  IIB  in  Figure  27 

f.  Case  IIC  in  Figure  30 
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5.   The  components  of  the  optimal  trajectory,  v, 
a,    0,  and  0,  together  with  the  related  airplane  parameters, 
rate  of  descent  and  normal  acceleration,  are  presented  for: 

a.  Case  IA  in  Figure  14 

b.  Case  IB  in  Figure  17 

c.  Case  IC  in  Figure  20 

d.  Case  IIA  in  Figure  25 

e.  Case  IIB  in  Figure  28 

f.  Case  IIC  in  Figure  31 

In  all  cases,  the  results  were  within  established 
specification  limits  and  were  considered  realistic  for  an 
airplane  landing.   The  anticipatory  nature  of  the  optimal 
control  was  evident  in  the  results,  especially  in  Case  IIA, 
where  the  optimal  control  produced  a  finite  thrust  pertur- 
bation at  the  initial  time,  although  the  airplane  was  as- 
sumed to  be  in  an  ideal  flight  condition.   This  apparent 
lack  of  continuity,  along  with  the  discontinuity  encounter- 
ed in  the  optimal  control  of  Cases  IB,  IC,  IIB,  and  IIC  at 
t  =  0,  was  predictable,  since  it  was  assumed  that  the  con- 
trols were  instantaneously  available.   Development  of  a 
state  model  to  include  control  dynamics,  presented  in 
Chapter  III,  will  presumably  eliminate  this  apparent  dis- 
crepancy. 
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c.    PROBLEM  AREA 

The  importance  of  proper  selection  of  the  constant 
values  for  the  elements  of  the  diagonal  weighting  matrices 
of  the  performance  measure  is  evident;   namely,  no  unique 
values  exist  and  hence  proper  selection  is  based  on  pro- 
ducing an  admissible  control  and  trajectory  when  the  op- 
timization technique  is  applied.   In  this  investigation, 
where  many  specifications  had  to  be  simultaneously  satis- 
fied, this  selection  became  paramount,  especially  when  it 
became  apparent  that  several  components  of  the  state  tra- 
jectory were  sensitive  to  changes  in  one  element  of  the 
weighting  matrices.   As  a  result,  the  trial-and-error  pro- 
cess required  to  obtain  values  for  these  weighting  matri- 
ces was  tedious  and  time-consuming.   Further  study  is 
warranted  toward  the  development  of  a  scheme  or  method  to 
aid  the  control  engineer  in  the  selection  of  these  weight- 
ing matrices. 
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VIII.   CONCLUSIONS 

For  the  results  presented  in  Chapter  VII,  the  follow- 
ing conclusions  are  stated: 

1.  The  feasibility  of  obtaining  realistic  results  by 
formulating  the  all-weather  landing  problem  as  a  linear 
tracking  problem  and  applying  optimal  control  techniques 
was  demonstrated. 

2.  Satisfactory  landings,  using  a  simple  mathemati- 
cal model  of  the  F-4J  airplane,  were  accomplished  in  all 
cases  investigated  when  the  derived  optimal  control  was 
applied. 

3.  Obtaining  values  for  the  elements  of  the  weighting 
matrices,  H,  Q,  and  R,  to  produce  admissible  and  realistic 
results  was  difficult  because  the  state  trajectory  was 
sensitive  to  changes  in  these  matrices.   The  process  was 
problematical,  since  a  large  number  of  design  specifications 
had  to  be  satisfied. 

Recommendations  for  logical  extensions  of  this  investi- 
gation follow: 

1.  Formulate  the  problem  to  include  the  control  dy- 
namics. 

2.  Formulate  the  problem  to  include  measurement  noise 
and  wind  effects. 

3.  Investigate  application  of  sub-optimal  control 
techniques  to  the  problem. 
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4.   Investigate  the  possibility  of  using  the  sensi- 
tivity of  the  state  trajectory  to  changes  in  the  weighting 
matrices  as  a  basis  for  developing  a  method  of  obtaining 
representative  values  for  the  elements  of  these  matrices. 
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APPENDIX  A    REPRESENTATIVE  DATA  FOR  THE  F-4J  AIRPLANE 
IN  THE  LANDING  CONFIGURATION 


Gross  Weight 

Equilibrium  Velocity  (v  ) 
StaU  Velocity  (V-uu) 
Equilibrium  Angle  of  Attack  (Ot  ) 


Stall  Angle  of  Attack  (£    ,, 

Limit  Sink  Rate  at  Touchdown 

Approximate  Thrust  Available 

MIL  engine  operation 

A/B  engine  operation 

Aerodynamic  Ground  Effect 

Control  Time  Constants 

Elevator  (Tg  ) 
e 


Thrust 


<V 


Dimensional  Stability  Derivatives 

X  =-0.593xl0_1sec"1         5 
v 

X  =  0.0 

q 

X  =  0.107xl02ft  sec2 
OL 

X&=  0.0 


Z6, 


32,500  lb 
223  ft/sec 
195  ft/sec 
0.33  rad 
0.49  rad 
23  ft/sec 

5,000  lb 
12,000  lb 
Negligible 

0.10  sec 
0.15  sec 


=  0.0 
=-0.132xl02ft  sec-2 


Z6T= 

Mv  = 


0.289xl0'3ft  lb"1sec"2 
0.450xl0"3ft"1sec"1 


-2 


x6=  0.0 
°e 

Xfi  =0.950xl0~3ft    lb"1sec' 

-2 

M     =-0.430    sec 

q 

Ma  =-0.130x10    sec 

Z     =-0.259    sec 

M^  =-0.260    sec 

Z     =0.0 

q 

Z    =-0.827xl02ft   sec"2 
Of 

_2 
M$  =-0.279  sec 

MA   =   0.0 
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APPENDIX  B    SUMMARY  OF  NUMERICAL  VALUES 


STATE  MODEL  CONSTANT 
a    =  -0.59300X10"1 


a12  =  0.10700x10' 
a13  =  -°-32172x10' 
a21  =  -0.11628x10 

a22  =  "0.37085 

a41  =   °-75232xl° 
a42  =  -0.12036x10 


-2 


-3 


a52  =  -°-22300x10' 


a__  =  0.22300x10 
b12  =  0.94667x10 
b21  =  -0.59193x10 
b22  =  0.12979x10 
b41  =  -0.27746x10 
b42  =  -0.33745x10 


-3 


-1 


-5 


-6 


a44  =  -0.69000 


c51  =  -0.11676x10' 


DIAGONAL  WEIGHTING  MATRIX  CONSTANTS 
CASE  I 


h,,  =  1.0 

h22  =  1.0 
h33  =  1.0 

h44  =  2«0xl° 
R    =  5.0 


-3 


qll  =  1*0xl° 

q22  =  1*0xl° 
q33  =  1.0x10 

q44  =  5.0x10 


-1 


-1 


-1 


-4 


CASE  II 


hll  =  5*0xl° 

h22  =  5*0xl° 

h  -  =  5.0x10 

h44  -  1.0 

hcc    =  5.0x10 

r   =  5.0 


-5 


-1 


-1 


-3 


qll  =  1-0xl° 


-5 


q22 

= 

l.OxlO"1 

q33 

= 

l.OxlO-1 

q44 

= 

5.0xl0_1 

q55 

= 

5.0xl0"4 

r22 

= 

5.0xl0"10 
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C  COMPUTER  PROGRAM 

C 

C 

C 

C  DEPTVATICN  OF  THE  OPTIMAL  CONTROL 

C         FOP  AN  ALL-WEATHER  AIRPLANE  LANDING  SYSTEM 

C 

C  NPGS  MASTER  THESIS,  JUNE  1969 

C 

C 

C  MAIN  PROGRAM 

C 

C 

REAL  *8  ITITLE( 12) »KKDGT(20) ,KK(20> ,XXDOT(5)  ,XX(5), 
♦  ITITLK36)  ,nTK,TOO,TFF,OTX 

REAL  *4  A(25),B(12),H(5),Q(5) ,R (4 ) , RRTF < 5  )  ,  AT ( 25  ) , 
*K(2C)  ,S(5),RRT(5)  .TEMPI  (25)  .TEMP2(25).TEMP3(25), 
*QC(2  5),SDCT(5) ,KDCT( 25 ) .X ( 5 ) , G( 5 ) .F(15).CC(5).P(15), 
*SVALS(5,1001),XVALS(5f1001),USTAR(2,100l),TBASE<1001), 
•  *BT(12),RIU),LAB(30),HH(5),HHH(5),TEMP4(25),Ll,L2tL3, 
♦GBASE < 600 ),ALFDCT( 1001) ,ALTDOT( 1001 ), ACCEL ( 1001 ) , 
*KVALS(15,10G1),XXX(1001)  ,  PSCAH600) 
C 

C     GRAPHICAL  OUTPUT  DATA. 
C       ITITLKI)  ARE  GRAPH  TITLES 
C       LAB(I)  ARE  CURVE  LABELS 
C 

RfcAD(5,524)  ( I T I TLl ( I ) , I =1 , 36 ) 

WFITE(6,5C5) 

WRITE  (6,  525)  (ITITLKI  ),I=1.36) 

READ(5,531)  (  LAB (  I  ) tI =1 , 15) 

READ(5,531)  (LAB ( I ), 1=16 .30 1 

WRITE* 6,507) (LAB(I),I=i,15) 

HRITE(6,5C7)(LAB(I ),I=16,30) 
C 

C     INITIALIZING  DATA 

C       NCASE  IS  CASE  NUMBER ( 1  =  1 .2  =  1 1 ) 
C       NO  IS  NUMBER  OF  PLANT  STATES 
C       NC  IS  NUMBER  OF  CONTROLS 
C       NSTK  IS  NUMBER  OF  INTERGRATION  STEPS 
C       DT  IS  INTEGRATION  STEP  SIZE 
C       TO  IS  INITIAL  TIME 
C       TF  IS  FINAL  TIME 
C       Ll,L2,L3  ARE  CONSTANTS  OF  OESIREC  TRAJECTORY  EGUATION 

READ(5,5CC)  NCASE .NG.NC. NSTK 

WPITE<6,501)  NCASE, NO, NC, NSTK 

PfcAD(5,5C2)  DT,T0,TF,L1,L2,L3 

DTX=DT 

DTK=-PTX 

WPITE(6,5C3)  DT,T0,TF,L1,L2,L3 

EPS=0.1E-0* 

NT  =  C 

NK=  NG*(N0+l)/2 

NKT=NK+NC 

N0N0=  NC*NC 

NCNC=  NG*NC 

NSTK1=NSTK+1 

NSTK2  =NSTK+2 

NSTK3=NSTK/2 

NSTK4=NSTK3+1 

NSTK5=NSTK+3 

RAD=57.295P 

GPAV=32.1725 

GAMMA=  -3.0/RAD 

VC=  223.0 

ALTD=100.C 
C 

C     TBASE  IS  TIME  BASE  FOR  COMPUTATION  WHILE  GBASE  IS 
C     TIME  BASE  FOR  GRAPHICAL  OUTPUT 


C 


TBASE(1)=0.0 
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DC  100  I=1,NSTK 
ICO  TEASE( I*1)=TBASE( II+DT 
DC  101  J=1,NSTK4 

101  GPASE(J)=  -TBASE(2*J-1) 
C 

C     PRGBLEM  DATA 

C       Ad  I  ARE  ELEMENTS  CF  A  MATRIX 

C       B(I)  ARE  ELEMENTS  CF  B  MATRIX 

C 

READ(5,5C4)  ( A( I )  ,1  =  1 ,25 > 

READ(5,504)  ( B ( I )  ,  1  =  1  ,10) 

DC  102  1=1, NO 

102  CCU  )=  0.0 
CC(NO)=VC*GAMMA 
IF(NCASE.EC2)  GO  TO  105 
DC  103  1=1, NO 

B(I)=  Bf 1*11 
A(I)=A(I+6) 

103  A(I+4)=  A(  I+il) 
DC  104  1=1, NO 
A(I48)=A(I+16) 

104  A(I+12)=A(H-21) 

105  WRITE(6,508) 
DO  106  J=1,N0 

106  WRITE(6,509)(A(I) , I=J,NCN0,N0) 
WRITE(6,51C) 

DO  107  J=1,N0 
10  7  WRITE (6,509) (B( I) , I=J,NCNC, NO ) 

WRITE(6,511) 

WFITE(6,516)   (CC(I).I=1,N0) 

PEAD(5,504)  (X(I)  ,1  =  1, NO) 

WRITE(6,523) 

WPITE(6,516)  (X(I),I=1,N0) 

READ(5,504)  (H( I  )  ,1  =  1  ,N0 ) 

REAC(5,504)  ( Q( I ) , 1=1 ,NG) 

READ(5,504)  (R ( I )  .  1  =  1  ,NC ) 

READ(5,504)  ( RRTF ( I ) , 1  =  1  ,N0) 

WPITE(6,512) 

WRITE(6,509)  ( H( I  )  ,1  =  1 ,NC ) 

WRITE(6,509)  ( 0( I ) , 1=1 ,NC) 

WRITE(6,5C9)  (R(I),I=1,NC) 

WRITE(6,522) 

WRITE(6,516)  (RRTF(I),I=1,N0) 
C 

C     FORM  A  TPANSPOSE,  B  TRANSPOSE,  AND  R  INVERSE 
C 


CALL  MTRA  ( A , AT ,NC  ,NC  ,0 ) 
CALL  MTPA(B,BT,NO,NC.O) 


IF(NC.EO.l)  GO  TO  109 
CALL  MSTR(P, TEMPI, NC, 2,1) 
CALL  SINV  (TEMPI. NC, EPS, IER) 
CALL  MSTR(TEMP1.RI,NC,1,2) 
IF(IER-l)  110,108,108 
103  WRITE(6,514) 
GC  TC  300 

109  RI(1)=    1.0/Rdl 

110  WRITE(6,526> 
DO    111    J=1,N0 

111  WRITE (6, 509) (AT (I ) , I=J ,NONO,NO) 
WRITE(6,527) 

DO    112    J=1,NC 

112  WRITE(6,509)(BT(I),I=J,N0NC,NC) 
WRITE(6,528) 
WRITE(6,509)(RI(I),I=1,NC) 

C 

C  FORM    FINAL    TIME    VALUES    FCR    K    AND    S 

C 

DO    113    1=1, NO 

113  RRTFU)=-RRTF(I) 
CALL    MSTR(H,K.NC,2,1) 

CALL    MPPD(H.RRTF,S,N0,NG,2,0,1) 
WRITE(6,515l 
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WRITE<6,516)  <K< I  )  ,  1  =  1, NK) 
WRITE<6,517) 

WPITEI6.516)  <  SU)  ,I  =  1,N0) 
DO  114  1=1, NO 

114  RPTU)=-PRTF(I) 
C 

C     FORM  P  MATRIX  WHICH  EGUALS  B*RI*BT 
C 

CALL  MPRD(B,RI .TEMPi.NC ,NC,0, 2 , NC ) 

CALL  MPPD<TEMPl,BT,TEMP2.NO,NC,0,0,NO) 

CALL  MSTR(TEMP2,P,NO,0,li 

WRITE(6,513) 

WPITE(6,516)<P(I)  ,1*1. NK) 

CALL  MSTR(0,QQ,NO,2,0> 
C 

C     FORM  KDOT  AND  SDCT  AND  INTEGRATE  BACKWARDS 
C 

DO  115  1=1, NK 

115  KK(I  )=K(I) 

DO  116  1=1, NO 

116  KMNK+I  )=S(I  ) 

DC  131   L=1,NSTK 
C 

C     THE  DESIRED  TRAJECTORY  IS  INTRODUCED  HERE.   IF  TO, 
C     TF,  DT.NSTK  OR  DESIRED  TRAJECORY  ARE  CHANGEC,  THE 
C     FOLLOWING  9  CARDS  MUST  BE  CHANGEC  ACCORDINGLY 
C 

IFCL.GE.401)  GC  TG  117 

DELT=  TBASE(401)-TBASE(L) 

TEMP=  EXP(L1*DELT)-1.0 

RPT<NO)=  (1.0/L1)*L2*TEMF  -  L2*DELT 

PRT(NO-l)=L3*DELT 

RPT(NO-2)=(L3/2.0)*DELT**2 

RRT(NO-3)=RRT(NC-2)-(L2/VO)*TEMP 

IF(NCASE.EO.l)  GO  TO  119 

RRT(N0-4)=  0,0 

GC  TO  119 

117  DO  118  1=1, NO 

118  RPT( 11=0.0 
115  DC  120  1=1, NK 

120  KVALSU  ,L)=K(I  ) 
DO  121  1=1, NO 

121  SVALS  (I,L)=S(I) 

122  CALL  MPRD(K.A,TEMP1,NC,N0,1,0,NC) 
CALL  MPRD  <AT,K,TEMP2,NC,N0,0,1,N0) 
CALL  MPRD  (K,P,TEHP3,NC,N0,l,ljN0) 
CALL  MPRD(TEMP3,K,TEVP4,NC,N0,0,1,N0) 
DC  123  I=1,N0N0 

123  KDOT(I)=  -TEMPKI  )-TEMP2 ( I ) -QQ< I ) *TEMP4<  II 
CALL  MSTRCKDOT, TEMPI, NO, 0,1) 

00    124    1=1, NK 

124  KKDOT(  I)=    TEMPKI  ) 

CALL  MSUB  (AT, TEMP3.TEMP4.N0, NO, 0,0) 
CALL  MPRD(TEMP4,S, TEMPI, NO, NO, 0,0.1) 
CALL  MPRD(0,RRT,TEMP2,NC,NO,2,0,1) 
DO  125  1=1, NO 

125  SD0T(I)=-TEMPKI)*TEMP2(I) 
DO    126    1*1. NO 

126  KKDOT(NK+I)=    SDOTU) 
SS=RKLDEQ(NKTfKK,KKOOT,TFF,DTK,NT) 
IFCSS-1.)    127,128,128 

127  WRITE(6.521) 
GO  TO  300 

128  DO  129  1=1, NK 

129  K(I)=  KKU) 
DO  130  1=1, NO 

130  SU)=    KMNK+I) 
IFISS-1.)     127,122t131 

131  CONTINUE 

DO    132    I=1,NK 

132  KVALSU.NSTKU)    =    K(I) 
DO    133    1=1, NO 
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133  SVALS(  I,NSTK+1)=S(I) 
C 

C     WRITE  K  AND  S  VALUES 


C 


NI=1 

IF(NK-IO)  134,134,135 

134  NR=NK 

NTI  =  1 
NRI  =  G 
GC  TO  136 

135  NTI=NK/10 
NR  =  1C 
NFI=NK-NTI*10 

136  DO  138  M=1.NTI 
WRITE(6,5ie) 

DO  137  J=1,NSTK1 

137  WFITEC6.509)  < KVALS< I  ,NSTK2-J) , I  =  NI , NR ) 
NI=NI+10 

138  NP  =  NP«-10 

IF(NRI)  141.141,139 
135  WRITE<6,518) 

NPI  =  NRH-NTI*10 

DO  140  J=1.NSTK1 
140  WPITE(6,509)  ( KVALS( I ,NSTK2-J ) , I  =  NI , NRI ) 

IC=0 


10  READ  (5.10C0,ENC=999)  CARD 
1000  FORMAT  (20A4) 


141  WPITE(6,52C) 
DC  142  J  =   UNSTK1 

142  WRITE<6,505)  ( SVALSC I  ,NSTK2-J ) , 1  =  1 , NO 
C 

C     OBTAIN  OPTICAL  CCNTROL  AND  OPTIMAL  TRAJECTCRY 
C 

NT=0 

DO  143  1=1, NO 

143  XX(I )=X(I) 
XXDOTCNO-3)  =  0.0 

XXDOT(NO)=  VO*GAMMA-»-VC*(X(NO-2)-X<NO-3)) 
DC  144  1=1, NO 

144  HHU)*0.0 

DO  157  L=1,NSTK 
DO  145  1=1, NO 

145  XVALS(I,L)=X(I ) 
C 

C     ALTDOT  IS  THE  RATE  OF  DESCENT  AND  ALFDOT  IS  THE 

C     DERIVATIVE  OF  ALPHA 

C 

ALTDOT(L)=XXDOT(NOI 

ALFDCT(L)=XXD0T(NQ-3) 

147  DC  148  1=1, NK 

148  TEMPHI)=  KVALS(I,NSTK2-L) 
DC  149  1=1. NC 

149  TEMP2(I)=  SVALSCI .NSTK2-L) 

CALL  MPPD(PI1BT,TEMP3,NC,NC12,0,NC) 

CALL  MPRD  (TEMP3, TEMPI, F  ,NC, NO, 0.1 ,NC) 

CALL  MPRD  (TEMP3,TEMP2,G,NC,NC, 0,0,1) 
C 

C     HH(NO)  IS  EQUILIBRIUM  ALTITUDE  TRAJECTORY  AND 
C     HHH  IS  F*HH 
C 

HH(NC)=ALTO  -( VC*3. C/RA0)*TBASE (L ) 

CALL  MPPD<F,HH,HHH,NC.NC, 0,0,1) 

CALL  MPRD(F,X,TEMP3,NC, NO, 0,0,1) 

DO  150  1=1, NC 
15C  TEMP4(I)=  -TEMP3(I)-G(I)  ♦HHH(I) 

DO  151  1=1, NC 

151  USTAR(I.L)=  TEMP4(I) 

152  CALL  MPRD(A,X,TEMPl,NC.NO, 0,0,1 ) 
CALL  MPRD  (B. TEMP4,TEMP2, NO, NC, 0,0,1) 
DO  153  1  =  1, NO 

153  XXDOT(I)=  TEMPI  UH-  TEMP2 < I )*CC < I ) 
SS=RKLDEQ(NO,XX,XXDOT,TOO,DTX,NT) 
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IFfSS-1. I  154,155,155 

154  WPITE(61521) 
GC  TO  300 

155  DO  156  1=1. NO 

156  X( !)=  XX(  I  ) 
IF(SS-1.)  154,152,157 

157  CGNTINUF 

DC  158  1=1, NO 

158  XVALS  (I.NSTK+1)  =X(I» 
ALTDOT(NSTK«-1)»XXDOT<NO> 
ALFDOT(NSTK+l)=XX00T(NC-3) 

C 

C     WRITE  OPTIMAL  CONTROL  AND  OPTIMAL  TRAJECTORY 

C 

16C  WPITE<6,529> 

DO  161  J=1,NSTK 

161  WRITE(6,530)   TBASE < J ) , (USTAR ( I , J )  ,  I*  1,NC  ) 
C     ACCEL  IS  THE  NORMAL  ACCELERATION 

DO  162  J=1,NSTK1 

162  ACCEL< J)=V0*(XVALS(N0-1,J)-ALFDCT(J) )/GRAV  *l.O 
WPITE<6,519) 

DO  163  J=1,NSTK1 

163  WPITE(6,53C)   TBASE ( J ) , ( XVALSl I  ,  J  1 , 1=  1  ,NC >  ,  ALTDOT ( J  ) , 
*ACCEL( J) 

C 

C     GRAPHICAL  OUTPUT SUBROUTINE  CRAW,DEVELCPEO 

C     AT  NPGS,  IS  USED  TO  CCNSTRUCT  THE  FOLLOWING  PLCTS 

C     (SUBROUTINE  AVAILABLE  AT  NPGS  COMPUTER  FACILITY) 

C  1.   PLOT  OF  OPTIMAL  CONTROL 

C  2.   PLOT  OF  OPTIMAL  AND  DESIRED  ALTITUDE 

C  TRAJECTORIES 

C  3.   PLOT  OF  VELOCITY.  ALPHA,  THETA,THETA  DOT, 

C  NORMAL  ACCELERATION  AND  RATE  OF  CESCENT 

C  4.   SOLUTION  TO  RICATTI  EQUATION 

C  5.   PLOT  OF  S  VALUES 

C 

DC  1  1=1.12 

1  ITITLEC I )  =  ITITL1(  I) 
DC  7  1=1, NC 
IF(NC.FO.l)  GO  TO  4 
GO  T0(2,3)  ,1 

2  MCC  =  1 
GFS=1.0 
GC  TO  5 

3  MCD=3 
GPS=C.0CC1 

GO    TO    5 

4  MOD=C 
GRS=1.0 

5  DC    6    L=1,NSTK3 

6  XXX(t )=USTAR( I  ,2*L-1)*GRS 

7  CALL    DRAW(NSTK3,XXX,GBASE,M0D,0,LAB( I ),ITITLE, 
*0. 05, 1.0, 11, 4, 2, 2, 9, 11, I, LA ST) 

DC    8    1=7.12 

8  ITITLE(I)  =  ITITL1U*6) 
DO    15    1=1,2 

GO    TO    (9, 13), I 

9  M0D=1 
C 

C     IF  DESIRED  ALTITUDE  TRAJECTORY  IS  CHANGED,  NEXT 

C     SIX  CARDS  MUST  BE  CHANGEC  ACCORDINGLY 

C 

DC  10  L=l,600 

10  XXX(L)=  ALT0*V0*GAMMA*T8ASE(L) 
DO  11  L=6C1,NSTK1 
DELT=TBASE(L)-TBASE<601) 
TEMP=EXP(L1*DELT)-1.0 

11  XXX(l)=ALT0*V0*GAMMA*TBASE(L)+L2*TEMP/Ll-L2*DELT 
DO  12  L=1,NSTK4 

12  XXX(L)=  XXX(2*L-1) 
GO  TO  15 

13  MCD=3 
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DO  14  L=1,NSTK4 

14  XXX(L)=  XVALS(N0,2*L-1) 

15  CALL  DRAW(NSTK4, XXX, G8ASE, MOD, 0, LABU+2)  ,  ITITLE, 
*2C.0t 1.0,1 l*li 2, 2, 9, 11,1, LAST) 

DO  16  1=7,12 

16  ITITLE(I)=  ITITLKI  +  12) 
DO  27  1=1,6 

GO  TO  (17, 18, 18, 18, 23,25), I 

17  M0D=1 
JJ*0 

IFCNCASE.E0.il  GO  TO  19 
GRS=0.01 

GO  TO  21 

18  M0D*2 
GRS*1.0 
GO  TO  21 

19  DO  20  L=1,NSTK4 

20  XXX(L)=0.0 
JJ  =  1 

GO  TO  27 

21  DO  22  L=1,NSTK4 

22  XXX(L)=  XVALS( I-J J  ,2*L-1 ) *GRS 
GO  TO  27 

23  DC  24  L=1,NSTK4 

24  XXX(L)=  ACCEL(2*L-1)*0.1 
GC  TO  27 

25  MCD=3 

DC  26  L=1,NSTK4 

26  XXX(L)=  ALTDOT(2*L-1>*0.01 

27  CALL  DRAH(NSTK4,XXX,GBASE,M0D,0,LAB(I+4) , ITITLE, 
♦0.05, 1.0, 11, 4, 2, 2, 9, 11,1, LAST) 

DO  28  1=7,12 

28  ITITLE( I)=ITITL1( 1418) 
DC  29  I=1,NSTK4 

29  GBASE(I)=  -TBASE(NSTK5-2*I) 
DC  37  1=1, NK 

DC  30  J=1,NSTK4 

30  PSCAL(J)=  KVALS(I,2*J-1) 

CALL  SCALE(NSTK4.PSCAL,60.0,GRS) 

GO  TO  (31, 32, 32, 32, 32, 32, 32, 32, 32, 33, 32, 32, 32, 32, 34), I 

31  M0D=1 

GO  TC  3  5 

32  M0D=2 

GC  TO  35 

33  IF(NCASE.EO.l)  GO  TO  34 
GO  TO  35 

34  M0D=3 

35  DO  36  L=1.NSTK4 

36  XXX(L)=KVALS(I  ,2*L-1)*GRS 

37  CALL  DRAW(MSTK4, XXX, GBASE, MOD, 0  ,LAB(H-10) ,  ITITLE, 
*60. 0,1. 0,1 1,4, 2, 2, 9, 11,1, LAST) 

DO  38  1=7,12 

38  ITITLE(M=  ITITLKI+24) 
DO  46  1=1, NO 

DC  39  J=1,NSTK4 

39  PSCAL(J)=  SVALSU  ,2*J-1) 

CALL  SCALE(NSTK4.PSCAL,3.0,GRS) 
GC  TO  (40, 41, 41, 42, 43), I 

40  MC0=1 

GO  TO  44 

41  M0D=2 

GC  TO  4* 

42  IF(NCASF.EO.l)  GO  TO  43 
GO  TO  44 

43  M0D=3 

44  DO  *5  L=1,NSTK4 

45  XXX(L)=  SVALSt I,2*L-1)*GRS 

46  CALL  DPAW(NSTK4, XXX, GBASE.MOD ,0 , LAB ( 1*25) , ITITLE, 
♦3.0,1. 0,1 1,4, 2, 2, 9, 11,1, LAST) 

300  CONTINUE 

500  FORMAT(4I10) 

501  FORMAT(    /// ,T8 , 'NCASE'1 ,14, /,T8, • N0=» , 14,/ , T8, 
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«"NC=',I4,/fT8, »NSTK=' ,14) 

502  FORMAT!6E10.5) 

503  FORMAT (T8, »0T= • ,E 13.5 , / , T8. 'TO=» ,E13. 5,/ , T8, ' TF=»  , 
*E13.5,/,T8,'L1=',E13.5,/,T8,'L2=',E13.5,/,T8, 
*'TF=',E13.5> 

504  FORMAT15E10.5) 

505  FORMAT!  '1'.//) 

506  F0RMAT(8E13.5) 

507  FORMAT(2X,15(A4,2X)  ) 

508  FORMAT!///, T9, 'A  MATRIX*,/) 

509  FORMAT! ICE13.5) 

510  FORMAT!///, T8,'B  MATRIX',/) 

511  FOPMAT(///,T8,'C  MATRIX**/) 

512  FORMAT!///, T8. 'H,Q, AND  R  MATR ICES ( DI AG  STORAGE)',/) 

513  FORMAT  !///,T8,'P  MATRIX  (SYMETRIC  STORAGE)1*/) 

514  FORMAT!///, T8, 'ERROR  IN  SINV',/) 

515  FORMAT!///, T8,«K!TF)  (SYMETRIC  STORAGE)',/) 

516  FORMAT(E13.5) 

517  FORMAT!///, T8,'S!TF)  VECTOR',/) 

518  FORMAT  ( ///, T8 , 'R ICATTI  COEF',/) 

519  FORMAT   ( ///,T8 ,' STATE  SOLUTIONS',/) 

520  FORMAT  !///,T8,»  S  COEF',  /) 

521  FORMAT  ( ///,T8 , 'ERROR  IN  PKLDEQ',/) 

522  FORMAT!///, T8,'RR!TF)  VECTOR',/) 

523  FORMAT!///,  T8,  »X!C)  VECTOR', /) 

524  F0RMAT(6A8) 

525  FORMAT(2X,(A8) 

526  FORMAT! ///,T8,  'A  TRANSPOSE',/) 

527  FORMAT!///, T8, '6  TRANSPOSE',/) 

528  FORMAT!///, T8, »R  INVERSE(DIAG  STORAGE)',/) 

529  FORMAT!///, T8, 'CPTIMAL  CONTROL',/) 

530  FORM*T(F6.2,8E14.5) 

531  FORMATQ5AA) 
END 

C 

c 

C  FUNCTION  RKLDEQ 

C 

C 

C  PURPOSE 

C  THIS    ROUTINE    SOLVES    A    SYSTEM    OF    N    FIRST-CROER 

C  ORDINARY  DIFFERENTIAL  EQUATIONS  BY  THE  RUNGE-KUTTA- 

C  GILL  FOURTH-CRCER  METHCD 

C 

C  DESCRIPTICN  OF  PARAMETERS 

C  N  -  NUMBER  OF  FIRST-CRDER  EQUATIONS  IN  SYSTEM  TO 

C  BE  SOLVED. (O.LE.N.LE. 25) 

C  Y  -  NAME  CF  LINEAR  ARRAY  OF  LENGTH  AT  LEAST  N, 

C  IN  WHICH  SCLUTION  VALUES  WILL  BE  STORED.   THE 

C  CALLING  PROGRAM  SHCULD  SUPPLY  INITIAL  VALUES 

C  BEFCPE  FIRST  ENTRY. 

C  F  -  NAME  OP  LINEAR  ARRAY  OF  LENGTH  AT  LEAST  N, 

C  IN  WHICH  THE  DERIVATIVES,  COMPUTED  IN  USER'S 

C  PROGRAM  APE  STOREC. 

C  X  -  THE  INDEPENDENT  VARIABLE,  WHICH  IS  ACVANCED 

C  WITHIN  RKLDEQ 

C  H  -  THF  INCREMENT  FOR  X,   WHICH  MAY  eE  CHANGED  AT 

C  THE  END  OF  ANY  INTERVAL.  (WHEN  S=2.0) 

c  M  -  AN  INTEGER  WHICH  COUNTS  THE  NUMBER  OF  TIMES 

C  ENTRY  TO  RKLDEQ  HAS  BEEN  MADE  DURING  CURRENT 

C  INTERVAL.   IT  MUST  BE  INITALLY  SET  TC  ZERO  BY 

C  USER  PtFORE  FIRST  CALL  OF  RKLDEQ.   SUBSEQUENTLY 

C  IT  SHCULD  NOT  BE  CHANGED  BY  USER. 

C  S  -  A  SWITCH  TC  BE  TESTED  BY  USER  UPON  RETURN  FROM 

C  RKLDEQ. 

C 

C  IF  S  =  1.0*  THE  CALLING  PROGRAM  SHCULD  NOW 

C  COMPUTE  VALUES  OF  F,  USING  CURRENT  VALUES  OF 

C  X  AND  Y,  AND  RETURN  TO  RKLDEQ. 

C  IF  S  =  2.0,  THE  END  OF  PRESENT  INTERVAL  HAS 

C  BEEN  REACHEO.   USER  SHOULD  STORE  ANC/OR  OUT- 
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C  PUT  CURRENT  X   AND/CR  Y  AND  TEST  FOP  ENC  OF 

C  COMPUTATION, 

C 

C 

FUNCTION  PKLDEO  ( N , Y ,F ♦ X ,H, NT ) 

KFAL*P  Y,F,X,H,G,H1,H2,H3,H6 

DIMENSION  Y(l)  ,  F(l),  0(25) 

NT  =  NT  +1 

GC  TO  (1,2, 3, 4), NT 

1  HI  =  H 

H2  =  HI  *  .5D0 
H3  =  HI  *  2. DO 
H6  =  H1/6.C0 
DO  11  J  =1,N 
11  0<J)  =  O.DC 
A  =  .5DO 
X  =  X  ♦  HZ 
GO  TO  5 

2  A  =  .2928932188134525 
GO  TO  5 

3  A  =  1.7071C67811865475 
X  =  X  ♦  H2 

GO  TO  5 

4  DO  41  I  =  1,N 

41  Y(I)  =  Y(I)  ♦  H6  *  F(I)  -Q(I)/3.00 
NT  =  0 
RKLDEO  =2, 
GO  TO  f- 

5  DO  51  L  =  liN 

Y(L)  =  Y(L)  ♦  A  *(H  *  F(L)  -0(D) 
51  OIL)  =  H3  *  A  *F(L)  ♦  ( 1. D0-3.D0*A )  *Q(L) 
RKLDEO  =1. 

6  RETURN 
END 

C 

C 

C  SLPRGUTINE  SCALE 

C 

C 

SUBROUTINE  SCALE ( NPTS  ,PTS ,SCLE ,GRS ) 

DIMENSION  PTSI600) 

RANGE=3.0*SCLE*SCLE/4.0 

MAX=ABS(FTS(1) ) 

DC  11  1=2, NPTS 

TMAX=ABS(FTS(I ) ) 

IF(MAX-TMAX)10,11,11 
1C  MAX=ABS(PTS( I) ) 

11  CONTINUE 
IF(RANGE-MAX)  12,22,13 

12  KFLAG=1 
IF(MAX.LE.(RANGE*10.0) )  GO  TO  14 
IF(MAX.LE. (RANGE*100.0) )  GO  TO  15 
IF(MAX.LE.(RANGE*1000.0) )  GO  TO  16 
IF(MAX.LE.(RANGE*1C000.0) )  GO  TO  17 
IFIMAX.LE. (RANGE*100000.0))  GO  TO  18 
GO  TO  18 

13  KFLAG=2 
IF(MAX.GE.(RANGE*0.1) I  GO  TO  22 
IF(MAX.GE.(RANGE*0.01)I  GC  TO  14 
IF(MAX.GE.(RANGE*0.001))  GO  TO  15 
IF(MAX.GF.(RANGE*0.0001) )  GO  TO  16 
IF(MAX.GE.<RANGE*0.00001)»  GO  TO  17 
IF(MAX.GE.(RANGE*0. 000001)1  GO  TO  18 
GO  TO  18 

14  TGRS=10.C 
GO  TO  1° 

15  TGRS=100.0 
GO  TO  19 

16  TGRS=100C.O 
GO  TO  19 

17  TGRS=10000.0 
GO  TO  19 
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18  TGRS=10COOC.O 

19  GO  TO  (20, 211 ,KFLAG 
2C  GPS=  1.0/TGRS 

Gt  TO  2  3 

21  GRS=  TORS 
GO  TC  23 

22  GFS=1.0 

23  CONTINUF 
RETURN 
END 
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